package com.asa.hutils.teizheng;

import java.util.Arrays;

import com.asa.hutils.teizheng.structure.Matrix;


public class Eigenvalue {
	
	/**
	 * 特征值
	 * @param mat
	 * @return
	 */
	public static Matrix getEigenvalues(Matrix mat) {
		Matrix hessenberg = convertToHessenberg(mat);
		int rows = hessenberg.rows();
		int cols = hessenberg.columns();		
		Matrix B = new Matrix(rows, cols);
		hessenberg.copyTo(B);
		int index = 0;
		while(true) {			
//			System.out.println(index++);
			if (index>1000000) {
				throw new RuntimeException("迭代次数过多，应该是算不出来了，可能解是虚数");
			}
			Matrix Q = Eigenvalue.getIdentityMatrix(rows);			
			Matrix R = new Matrix(rows, cols);
			B.copyTo(R);
			for(int i=0; i<rows-1; i++) {
				Matrix partR = getRForAStep(R, i);
				R = partR.multiply(R);
				Matrix RTranspose = partR.t(); 				
				Q = Q.multiply(RTranspose);
			}	
			B = R.multiply(Q);
			if(isConvergence(B)) {
				double[] b = new double[cols];
				Arrays.fill(b, 0);
				double[] eigenvalue = getEigenvalue(B);	
				Matrix eigenvalueMat = new Matrix(eigenvalue.length, 1, eigenvalue);
				return eigenvalueMat;
			}
			index++;
		}
	}

	private static double[] getEigenvalue(Matrix B) {
		double[] eigenvalue = new double[B.rows()];
		for(int i=0; i<eigenvalue.length; i++) {
			eigenvalue[i] = B.at(i, i);
		}
		return eigenvalue;
	}

	private static Matrix getRForAStep(Matrix R, int i) {
		Matrix partR = Eigenvalue.getIdentityMatrix(R.rows());
		double r = Math.sqrt(R.at(i, i)*R.at(i, i) + R.at(i+1,i)*R.at(i+1,i));
		partR.set(i, i, R.at(i, i) / r);
		partR.set(i, i+1, R.at(i+1, i) / r);
		partR.set(i+1, i+1, R.at(i, i) / r);		
		partR.set(i+1, i, -1*R.at(i+1, i) / r);
		return partR;
	}
	
	private static Matrix convertToHessenberg(Matrix A) {
		int rows = A.rows();
		int cols = A.columns();
		Matrix R = new Matrix(rows, cols);
		A.copyTo(R);
		HessenbergHouseholderReflection hessenbergHouseholder = 
				new HessenbergHouseholderReflection();
		for(int j=0; j<rows - 2; j++) {
			Matrix H = 
				hessenbergHouseholder.getHouseholderReflectionAtACol(R, j);
			if(H.rows() == 0) {
				return R;
			}			
			R = H.multiply(R);
			R = R.multiply(H);		
		}		
		return R;
	}
	
	private static boolean isConvergence(Matrix B) {
		final double exp = 0.0001;
		for(int i=1; i<B.rows(); i++) {
			if(Math.abs(B.at(i,i-1)) >= exp) {
				return false;
			}
		}
		return true;
	}

	private static Matrix getIdentityMatrix(int size) {
		Matrix I = new Matrix(size, size);
		for(int i=0; i<I.rows(); i++) {
			for(int j=0; j<I.columns(); j++) {
				if(i == j) {
					I.set(i, j, 1);
				} else {
					I.set(i, j, 0);
				}
			}
		}
		return I;
	}
	
}
